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1.  Introduction 


Piezoelectric  materials  possess  the  unique  property  that  when  subjected  to  a  deformation,  they 
produce  electricity.  Conversely,  when  these  materials  are  subjected  to  an  electric  field,  they 
undergo  defonnation.  This  unique  property  has  been  exploited  by  researchers  and  engineers  to 
achieve  desired  effects  such  as  in  the  design  of  the  ignition  system  in  the  rocket-propelled 
grenade  launcher  (RPG)  (shown  in  figure  1);  it  continues  to  intrigue  researchers  with  its  promise 
in  the  design  of  new  systems  and  devices.  New  piezoelectric  applications  are  limited  only  by  the 
imagination.  Figure  2  shows  rotor  blade  research  conducted  at  the  University  of  Maryland  in  an 
attempt  to  design  smart  controllable  blades  that  have  desirable  twist  characteristics. 


Source:  http://science.howstuffworks.com/rpg3.htm. 
Figure  1.  RPG  launcher. 


Source:  http://www.enae.umd.edu/AGRC/ 
pict/. 

Figure  2.  Smart  controllable  rotorcraft  twist 
research. 
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At  the  United  States  Military  Academy,  students  and  researchers  (Labo,  2005)  have  researched 
the  incorporation  of  piezoelectric  materials  in  combat  boots  (see  figures  3  and  4).  This  research 
was  motivated  by  an  attempt  to  capture  the  excess  energy  generated  during  the  natural  gait  of  a 
soldier’s  step.  Other  uses  of  piezoelectric  materials  include  sophisticated  power  generation 
systems  for  microsensors  and  microelectromechanical  systems  devices  and  household  uses  such 
as  ignition  switches  for  propane  barbeque  grills. 


Figure  3.  Piezoelectric  board  inserted  into  a  combat  boot. 


Figure  4.  Piezoelectric  board. 


2.  PZT  Monoclinic  Structure 


The  introduction  to  piezoelectric  materials  just  mentioned  motivates  an  explanation  of  how 
piezoelectric  materials  work.  To  do  this,  the  molecular  makeup  of  these  materials  must  be 
examined.  To  limit  the  scope  of  this  explanation  and  stay  within  the  context  of  the  available 
research,  the  piezoelectric  material  PbZr052Ti04&( 93  will  be  discussed.  The  material 
PbZr0  52Ti04SO2  (see  figure  5)  is  a  subset  of  piezoelectric  materials  called  PZTs  that  contain 
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Figure  5.  The  PbZrQ52TiQ4&03  (PZT)  molecule  (Rabe,  2002). 

lead  (Pb),  zirconium  (Zr),  and  titanium  (Ti).  This  material  is  also  a  subset  of  materials  called 

perovskites  that  are  composed  of  atoms  of  the  ABO}  form.  In  PbZr052Ti04SO3,  Pb  is  the  A  atom 

and  either  Zr  or  Ti  (randomly  distributed)  is  the  B  atom. 

The  author  of  figure  2  (Rabe,  2002)  aptly  describes  this  molecule’s  behavior  as  follows: 

The  ideal  perovskite  structure  (ABO3)  centred  on  a,  the  B  site  (Ti/Zr) 
and  b,  the  A  site  (Pb).  In  the  piezoelectric  material  PZT  the  Zr  and  Ti 
atoms  are  randomly  distributed  on  the  B  sites.  Using  first-principles 
calculations  Grinberg  and  colleagues  show  that  each  Pb  atom 
displaces  preferentially  towards  square  green  faces  with  the  largest  Ti 
fraction.  These  local  distortions  can  be  correlated  with  global  structure 
and  behaviour.  For  example,  in  PZT  the  displacement  of  the  metal  ions 
(Pb,  Ti  and  Zr)  away  from  the  centre  of  their  oxygen  cages  generates 
local  polarization  and  ultimately  ferroelectricity.  c,  The  net 
polarization  and  overall  structure  of  PZT  changes  with  the  ratio  of 
Ti:Zr  atoms,  as  shown  in  th ex~T phase  diagram.  The  six  structural 
phases  of  PZT  at  ambient  pressure  are  mostly  ferroelectric,  with  the 
exception  of  antiferro electric  behaviour  (A)  near  PbZrCU  and  a 
paraelectric  cubic  phase  (C)  at  high  temperatures.  Grinberg  et  al.  show 
that  the  complex  phase  behaviour  of  different  Zr/Ti  compositions  can 
be  explained  by  averaging  the  changing  distribution  of  Pb 
displacements,  which  depend  on  their  local  environment.  T  =  Ti-rich 
tetragonal  phase;  R  =  rhombohedral  phase;  M  =  monoclinic 
ferroelectric  phase. 
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3.  Notation 


The  following  portion  of  this  report  concerns  itself  with  incorporating  the  piezoelectric  behavior 
of  the  PZT  molecule  shown  in  figure  5  into  a  mathematical  fonnulation  which  can,  in  turn,  be 
implemented  by  computational  methods.  It  is  appropriate  at  this  point  to  make  clear  distinction 
of  the  notation  used  throughout  the  rest  of  this  report.  To  discern  between  scalar  and  tensoral 
variables,  we  will  use  the  standard  indicial  notation  ut  for  first-order  tensors  (i.e.,  vector 
variables),  u:j  for  second-order  tensors,  uikj  for  third-order  tensors,  etc.  In  a  three-dimensional 
(3-D)  calculation,  all  tensor  and  vector  indices  span  the  range  (1,2,  and  3).  Variables  with  no 

indices  will  be  scalar  quantities.  Moreover,  variables  in  bold  font  (e.g.,  u)  will  be  tensor 
variables  (order  implied  by  usage),  and  unbold  font  will  imply  scalar  quantities.  When 
discussing  atomistic  relationships,  atom  labels  will  be  denoted  by  the  lone  indices  i,  j,  or  k  in 
parentheses,  usually  as  superscripts.  For  example,  r(y)  (or  r (i,j)  )  is  the  position  vector  from 
atom  i  to  atom  j  (see  figure  6). 


Figure  6.  Three  atoms  depicted  in  the  reference  and  current  configurations. 
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In  this  report,  the  term  “fine  scale”  refers  to  events  occurring  at  atomistic  scales  on  the  order  of 
10  Angstroms,  while  the  tenn  “macroscale”  refers  to  events  occurring  at  scales  larger  than 
100  pm  (i.e.,  visible  to  the  naked  eye).  The  term  “coarse  scale”  has  a  meaning  synonymous  with 
macroscale. 


4.  Atomistic  Behavior  and  the  Constitutive  Equation 


In  fonnulating  mathematical  equations  to  model  the  physical  world,  material  properties  are 
incorporated  via  constitutive  equations.  The  classical  constitutive  equation  used  in  solid 
mechanics  on  the  macroscale  is  Hooke’s  law,  and  the  tensor  variable  possessing  the  material 
behavior  is  the  elasticity  tensor  denoted  CIJKL .  Hooke’s  law  is  stated  as  follows: 


S !J  C ' ijkl^kl  ’ 


(1) 


where  Su  is  the  stress  tensor  and  EKI  the  strain  tensor.  In  a  loose  sense,  Su  embodies  the 
forces  acting  on  a  body  while  E ^  embodies  the  displacements  that  occur  as  a  result  of  these 
forces.  It  is  a  reasonable  assumption,  therefore,  that  if  material  behavior  normally  associated 
with  the  PZT  fine  scale  is  to  be  mathematically  embedded  into  equations  describing  the 
macroscale,  the  elasticity  tensor  seen  in  Hooke’s  law  is  where  the  embedding  should  occur.  The 
formulation  of  an  elasticity  tensor  encapsulating  PZT  atomistic  behavior  is  needed  for  research. 


5.  Atomistic  Elasticity  Tensor 


5.1  Atomistic  Kinematics 


Figure  6  depicts  three  atoms  (i,  j,  and  k)  in  the  reference  and  current  configurations.  These 
atoms  are  located  by  Cartesian  position  vectors  Xa/),X(/A),X(/A)and  \(lJ> ,\<lk) ,  and  x(jk)m  the 
reference  and  current  configurations,  respectively.  The  atoms  may  locally  be  described  by  their 
proximity  to  each  other  by  the  position  vectors  R(//),R(/A),R</A)  and  r(iJ\  rhk> ,  and  r(jk> ,  also 
shown  in  respective  configurations. 

The  angle  created  by  two  vectors  subtending  atom  i  is  defined  by  the  symbol  0(jlk)  =  0a'J) ,  where 
the  atom  subtended  is  the  middle  index.  The  cosine  of  this  angle  is  therefore  as  follows: 


cos# 


urn  _ 


.(ij) 


.(*) 


.an 


M) 


y(‘j)  y(ik) 


(2) 


and 


r(ij)  —  .  y(’j) 


(3) 
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We  note  that  all  kinematic  variables  defined  for  any  three  atoms  can  be  simplified  to  two  vectors. 
For  example,  0(jlk)  =  0,j'k’(r<ij> ,rllk>)  and  rijk]  =  ri,k)  -  r(ij) .  Therefore,  in  what  follows,  the  two 
truly  independent  variables  are  r<iJ)  and  r(lk> . 

5.2  Definition  of  Potential 

The  fonn  of  the  potential  is  given  as  follows: 

<!>=  Hv2(i,j)+  Hv3(i,j,k),  (4) 

i,j  iJ,k 

i<  j  i<  j<k 


where  the  first  summation  denotes  pair  interactions  and  the  second  summation  the  interactions  of 
atom  triples.  Note  that  there  is  no  factor  of  1/2  contained  in  the  equation  because  the  summation 
conditions  disallow  double  counting.  The  pair  and  triple  terms  are  defined  as  follows: 

v2(i,j)  =  sf2(r(ij> /cr),  (5) 

and 

v2(i,j,k )  =  s  f3( r(iJ)  / <j,x"k>  /  cr,x(jk)  /  <j)  ,  (6) 

where  the  constant  s  is  chosen  to  give  the  minimum  of fi  at  -1,  and  cr is  chosen  to  make  /2(2  '' ) 
vanish.  These  constants  are  given  for  Si  in  (Stillinger  and  Weber,  1985). 


The  pair  term  is  a  function  of  the  scalar  distance  between  two  atoms,  i  and  j,  and  is  defined  as 
follows: 


| A(Br  p -r  9)exp[(r-a)  1  ], 


r  <a 

? 

r  >  a 


(7) 


where  A,  B,  p,  q,  and  a  are  constants.  The  range  of  the  pair  interactions  is  limited  to  a  circular 
region  of  radius  a  around  each  atom.  The  triple  interaction  is  defined  as  follows: 

/3(r®,rW,rW)  =  h(/iJ\rm,0Uik))  +  h(rui\ruk),0m)  +  h(rm ,rm ,0m) ,  (8) 


where 

h(rm,r(ik),0um)  =  /lexp [r(rUJ)  -  a)  '  +  y(r(lk)  -  a)  '  Jcos<9r  /'Vf)  +  .  (9) 

Based  on  the  previous  discussion  that  there  are  only  two  true  independent  variables,  we  can 
redefine  h  in  the  fonn  of  some  functions  h (,)  and  hU)  so  that  equation  8  can  be  rewritten  as 
follows: 

/3  =  h{,) (|r(,;?) |,|r(tt) I)  +  h  U) (|r°'° |,|r w I)  +  h{k) (|r(fa'} I, |rW) I) 

=  hl,)  (rm,rw)  +  hU)  (rao,rw)  +  h{k\ r(fe),rW)) 

=  h®(rmJm)  +  hU)(-rm,rm +  ,rm  -rw),  (10) 
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where 


hU)(r{ij)  ,ruk))  =  Texp 


y  *Jr(ij)  ■  r(iJ)  -  a)  +  y  Vr(,i)  •  rw 


-(»>■) 


.(y) 


.(*) 

+  1 

.(*)!  3 


(ID 


5.3  Finite  Strain  and  Classical  Continuum  Mechanics 

The  atoms  labeled  i,j,  and  k  are  assumed  in  the  defonned  configuration.  The  equations  thus  far 
are  only  considering  the  instantaneous  configuration  to  evaluate  the  energies  from  equation  4. 
From  a  Lagrangian  continuum  solid  mechanics  perspective,  it  is  important  to  know  how  those 
atoms  got  there.  We  therefore  define  atoms  labeled  /,  J,  and  K  as  the  reference  configuration 
relative  to  the  original  Cartesian  basis. 

We  consider  only  affine  deformations,  meaning  that  the  deformation  process  at  the  atomic  scale 
is  assumed  to  work  in  the  same  way  as  the  continuum  scale.  Under  this  one-to-one  mapping 
assumption,  position  vectors  in  the  reference  configuration  R  get  transfonned  to  the  current 
configuration  r  through  the  deformation  tensor  F, 

(r);  =  (FR),  =  r,  =  J^FuR,  =  F^R,  ,  (12) 

/ 

where  repeated  indices  are  assumed  summation. 

We  use  the  conventional  strain  measures  for  the  Green  strain  and  Right  Cauchy-Green  strain, 
respectively,  as  follows: 

EiJ=\(F,iFu-8u), 


and 


Cu=FiIFu.  (14) 

We  can  then  obtain  a  host  of  elasticity  tensors  by  taking  second  derivatives  of  a  strain  energy 
function.  From  equation  4,  the  strain  energy  function  is  defined  by  the  following: 

(15) 


where  N  is  the  number  of  atoms.  Different  types  of  elasticity  tensors  may  then  be  detennined. 
For  example,  we  may  choose  among  the  following  fourth-order  tensors: 


a2o  a2o  a2o  a2o 

dCndC^  ’  dEudEKL  ’  dFudFkL  ’  dsijdskl  ’ 


(16) 
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where  the  last  tensor  is  based  on  the  infinitesimal  strain  and  would  give  the  standard 
infinitesimal  elastic  constants  after  taking  appropriate  directional  derivatives. 


The  process  for  obtaining  any  of  the  elastic  property  tensors  is  identical.  We  make  repeated  use 
of  the  chain  rule  to  analytically  detennine  the  derivatives.  Let  us  consider  the  following  elastic 
tensor  as  the  illustrative  example: 


' IJKL 


d20 

dCu8Ca 


(17) 


By  making  use  of  earlier  definitions,  we  proceed  as  follows: 

820 


" IJKL 


dCudCn 

82  O 


8FmN  dFpQ 


8FmN8FpQ  8CU  8C 


<92®  8r,(iJ)  8r 


m 


-+- 


82 ®  8r,(iJ)  8r 


(y) 


8r™dr™  8FmN  8F pQ  r/f'r/f  8FmN  8F pQ 


8FmN  8F, 


pQ 


+ 


32®  8rm  8r?k)  a2®  8r(ik)  8rf 


8r™8r™  8FmN  8F pQ 


+ 


8CU  8C a 


8rw8r^ 


8FmN  8FpQ  j 


(18) 


It  is  important  to  recall  that  the  summation  convention  does  not  apply  to  indices  in  parentheses. 
To  complete  the  process  of  equation  18,  we  use  the  definition  in  equations  13  and  14  to  obtain 
the  following: 


8F, 


pQ 


8C 


KL 


(scF 

-1 

{ 

i  8F*>  J 

-1 


pL^QK  +  F pK 8 Q,  t 


(19) 


and  equation  12  to  obtain 


dr™  _8(FiIR^))_ 


8F„, 


8F„, 


S  8  RUJ>  =  8  R 

UmiUNI1VI  Umi1VN 


00 


(20) 


The  only  task  that  now  remains  is  to  evaluate  the  derivatives  of  the  strain  energy  potential  ® 
with  respect  to  the  atomic  vectors.  This  is  also  a  straightforward  procedure.  We  start  with  the 
first  derivative  as  follows: 


f 


d</> 

dr(p 


d 

8rp^ 


Xv2  (*’■)) + 

ij  Uj,k 

\ z< j  i< j<k 


8 


(21) 


Z_^2_+  y  ^v3 
a  (i/)  Zj  a  (y) 

<J  0VP  U,k  °rp 

i<  j  i<  j<k 

df2  0r&) 


£ 


=  ?  £ 


£ 


^  0/j(O  d/l00  0/?W  A 


i,j  5r<V)  5rpy)  ij.k 
i<  j  i<  j<k 


ydrpJ)  dr ^  dr^  j 


where  the  derivative  in  the  first  term  is  given  by  the  following: 


dr(ij)  dJr  r 


.(y)„(y) 


1 


dr(y  d/u) 


=  — L  •  2r}ij)S,„  = 


Av) 


‘p  rm  • 


(22) 


Note  that  drUj)  /  dr‘k]  =  /  dr  'J)  =  0  .  Also,  the  derivative  of  the  nonzero  part  of/)  is  as 

follows: 


f^  =  ^exp[(r-a)-1] 
or 


-Br 


-p 


P+- 


1 


r  ( a  -  r) 


+  r 


U- 


1 


r  ( a  -  r)  ) 


(23) 


A  similar  derivative  needs  to  be  evaluated  for  hU) . 


Likewise,  the  key  second  derivatives  that  need  to  be  derived  are  as  follows: 

d2(f)  d2(j)  d2</>  d2</> 

dr^dr^  ’  dr^j)d/qik)  ’  dr(‘k) 7drf)  ’  dr^dr™  ’ 


which  will  contain  expressions  such  as 

3rm  _d(r™/r™)_  1  M)d(l  lrm) 

QrW)Qr(‘j)  ~  Qr(ij)  ~  r(v)  °pq  +  rp  Qrm 


r 


I A  01  ..'01  N 

+  _P _ 

}A  01  jAOl 


(24) 


(25) 


By  observing  the  equations  thus  far,  the  elastic  stiffness  tensor  is  fully  determined.  The 
remaining  task  is  extensive  symbolic  manipulations  to  obtain  a  closed  form  expression  (if  so 
desired)  of  the  tensor,  with  repeated  use  of  partial  derivatives  and  chain  rule  multiplications. 
Due  to  the  algebraic  complexity  of  the  potential  function,  the  final  expression  is  generally 
lengthy,  but  still  obtainable. 


6.  Variational  Multiscale  Formulation 


While  the  elasticity  tensor  formulated  in  section  5  will  confine  the  displacement  of  PZT  atoms  to 
proper  paths  and  distances  within  the  atomic  lattice,  a  coupling  between  these  atomic 
displacements  and  the  electric  field  E  is  needed.  Moreover,  the  governing  equations  for  E 
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(i.e.,  Maxwell’s  equations)  must  be  present  in  the  overarching  mathematical  model  since  they 
share  a  role  equal  to  those  of  solid  mechanics  regarding  piezoelectric  phenomenon.  To  allow  for 
this,  the  equations  used  in  this  model  will  be  formulated  in  a  variational  multiscale  framework. 

6.1  Multiscale  Deformation 

In  the  variational  multiscale  context,  the  displacement  u  of  any  point  in  the  body  (assumed 
synonymous  with  an  atom’s  displacement)  is  decomposed  into  two  parts  u“  and  u  r\  whereby 
u“  is  a  local  displacement  on  the  fine-scale  and  ifi  is  a  global  displacement  on  the  macroscale. 
Figure  7  shows  a  pictorial  example  of  this  notion. 


Figure  8  demonstrates  some  additional  key  points  regarding  variational  multiscale  defonnation. 
In  this  figure,  a  hypothetically  thin  rod  is  used  to  demonstrate  that  (1)  displacement  u"  will  have 
a  magnitude  on  the  order  of  Angstroms,  while  u77  will  have  a  magnitude  on  the  order  of  meters; 

(2)  strain  is  scale  independent  -  a  point  that  will  be  very  important  in  subsequent  sections;  and 

(3)  displacement  u/;  occurs  due  to  the  summation  of  a  very  large  number  of  fine-scale 
displacements  u“  occurring  over  the  entire  body. 
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Thin  rod 
in  tension 


1 

0.01  meter 


1  Angstrom 


T 


1  meter 


Atoms 

(not  to  scale) 


F 


KEY  POINT  #1:  The  magnitude  of 
il°  is  on  the  Angstrom  scale.  The 
magnitude  of  UP  is  on  the  meter 
scale. 


Local  strain  is  0.01 


Angstrom 

Angstrom 


KEY  POINT  #2:  Strain  is  scale 
independent.  It's  equal  to  0.01 


Global  strain  is  0.01  -^Hr 


KEY  POINT  #3:  An  enormous 
number  of  u°  experienced  by  all 
the  atoms  over  a  region  sums  up  to  up  . 


Figure  8.  Some  key  points  made  by  observing  the  deformations  in  a  hypothetical  thin  wire. 


6.2  Multiscale  Kinematics 

To  accommodate  fine-scale  and  coarse-scale  equations,  an  additional  intermediate  configuration 
P  is  introduced  to  the  two-configuration  arrangement  seen  in  classical  continuum  mechanics  (see 
figure  6).  The  three  configurations  are  shown  in  figure  9.  In  this  figure,  the  position  of  a  point 
in  the  intermediate  configuration  is  denoted  by  X  . 


With  these  configurations,  a  logical  definition  of  variables  leads  to  simplistic  equations.  If 
displacements  are  defined  as 

u“  :=  X -  X  ,  (26) 

up  :=  x - X,  (27) 

and 


then,  as  desired,  u  =  u“ 


u  :=  x  -X, 

+  up .  Defining  the  deformation  gradients  as 

F:=^, 

dX 


(28) 


(29) 
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Figure  9.  Multiscale  kinematics.  Three  configurations  are  used  in  this  framework  - 
reference,  intermediate,  and  current  (Creighton  et  al.,  2004). 


and 


yields  the  following  equations: 


and 


F  =  FPFU, 

du° 


F“  =  I  +  ■ 


0X 


Fp  =  I  -  Vup , 


(30) 

(31) 

(32) 

(33) 

(34) 


where  the  nabla  symbol 
d(*) 

such  that  V(*)  := - . 

dx 


represents  derivatives  taken  with  respect  the  current  position  vector 
From  these  equations,  we  define  the  multiscale  Green  tensor  as 


x 


Ep  :=-(Fp,Fp  +1 

2  v 


(35) 


and  the  multiscale  Right  Cauchy-Green  tensor  as 
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CP  J7PTJ7P 


(36) 


For  a  rigorous  derivation  of  the  variables,  see  Creighton  et  al.  (2004). 


6.3  Multiscale  Governing  Equations  of  Solid  Mechanics 

With  the  kinematic  variables  defined,  we  recall  macroscale  equations.  Equations  37-39  are  (1) 
the  conservation  of  linear  momentum,  (2)  a  push  forward  of  the  second  Piola-Kirchhoff  stress 
tensor  S  to  the  current  configuration,  and  (3)  Hooke’s  law.  Note  the  appropriateness  of  the  J3 

superscript  in  the  following  macroscale  equations  (recall  that  up  is  macroscale  displacement): 

V  •  tr  +  b  =  piip ,  (37) 

o  =  -  FPSFP\  (38) 

J 

and 

S  =  C  :  Ep .  (39) 


Recalling  equation  17  from  section  5,  the  multiscale  version  of  this  equation  becomes 

d2<s> 

~  dCpdCp  ' 


(40) 


It  should  be  noted  that  using  C/;  in  the  derivation  of  the  atomistic  elasticity  tensor  C  is 
appropriate  since  strains  are  scale  independent  (see  key  point  no.  2  in  figure  8).  Conversely, 
incorporating  the  macroscale  displacement  variable  up  into  the  formulation  of  C  would  be 
inappropriate  due  to  its  scale. 


6.4  Governing  Equations  of  Electromagnetism 

The  governing  equations  for  electromagnetism  are  Maxwell’s  equations.  From  these  equations, 
a  relationship  for  the  electric  field  vector  E  can  be  written  as  follows: 

—  V  x  V  x  E  +  eE  +  yE  =  0  .  (41) 

(See  Weidlinger  Associates  Inc.  [1991]  for  the  derivation.)  Recognizing  that  the  force  exerted 
on  a  particle  (or  atom)  of  charge  q0  in  an  electric  field  is 


F  =  qo  E  (42) 

and  that  Newton’s  Second  Law  for  the  atom  is 

F  =  mi“,  (43) 
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these  two  equations  can  be  combined  to  create  the  following: 


mu 


a 


or,  alternatively, 


where 


m 


and  m  = 


(0.48  •  mass(Ti )  +  0.52  •  mass(Zr)) 

2 


for  the  PZT  at  hand. 


(44) 


(45) 


Consequently,  equations  37,  41,  and  45  constitute  three  equations  and  three  unknowns  which 
adequately  describe  the  problem. 


6.5  Variational  Equations 

Casting  equations  37,  41,  and  45  into  the  weak  fonn  for  a  numerical  solution  by  the  finite- 
element  method,  the  problem  is  expressed  mathematically  as  follows: 


:=  J  Vwp :  o(uu,  u|!)dSt— J  (wp  •  b)dSt-  J (wp  •  t)d5St-J  (wp  •  pup)dSt  =  0  (46) 


as, 


with  unknowns  {ua,up}. 


^  :=  J  wa  •  (/ui“  -  qE)dSt  =  0 


(47) 


with  unknowns  {ua ,  E}  . 


^3  :=  J  w" 


A 


—  V  x  V  x  E 

H 


dSt  +  J  wE  •  (eE  +  yE)dSt  =  0 


(48) 


/  s, 


with  unknown  {E} . 


7.  Conclusions  and  Future  Work 


7.1  Atomistic  Elasticity  Tensor 

The  elastic  property  tensor  will  be  a  3  x  3  x  3  x  3  fourth-order  tensor  which  will  have  8 1 
nontrivial  tenns.  All  trivial/nontrivial  terms  and  symmetries  associated  with  the  elastic  tensor 
are  determined  solely  by  the  crystal  and  NOT  prescribed  a  priori.  The  challenge  for  the  PZT 
structure  at  hand  will  be  to  capture  the  behavior  of  the  different  types  of  atoms  (which  vary 
significantly)  when  formulating  the  elasticity  tensor  C  .  Moreover,  using  the  potential  in 
equation  9  is  meant  merely  as  an  illustrative  example  of  the  procedure  since  it  neglects  long- 
range  interactions  common  in  piezoelectric  materials. 
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7.2  Multiscale  Formulation 

Although  the  equations  are  all  valid,  issues  concerning  stability  have  not  been  addressed  and  will 
not  be  apparent  until  the  stated  weak  fonn  equations  have  been  implemented  and  tested. 
Boundary  conditions  for  u“  and  up  also  need  to  be  addressed  but  can  be  modeled  from  earlier 
work  in  the  field.  An  initial  approach  might  be  to  use  periodic  boundary  conditions  for  u“  (used 
in  most  fine-scale  molecular  dynamics  simulations),  and  boundary  conditions  typical  of  classic 
macro-scale  solid  mechanics  for  up . 
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